##############################################################
#   S-puls program (R program) for computing the confidence   
# coefficient of the 95% Agresti-Coull interval. 
#   n is the sample size.                    
##############################################################   

n <- 5
k <- 1.96
n1 <- n+k^2
x <- seq(0,n,by=1)
x1 <- x+k^2/2
p <- x/n
q <- 1-p
p1 <- x1/n1
q1 <- 1-p1
ci <- cbind(p1-k*sqrt(p1*q1)*n1^(-1/2),
            p1+k*sqrt(p1*q1)*n1^(-1/2))

aa <- matrix(sort(ci))
bb <- matrix(aa[aa>0])
cp <- matrix(bb[bb<1])


datalist <- array(0,c(nrow(ci),ncol(ci),nrow(cp)))
for(k in 1:nrow(cp))
{
 datalist[ , ,k] <- cbind(cp[k,1]-ci[ ,1],ci[ ,2]-cp[k,1])
}

index <- matrix(0,nrow(ci),nrow(cp))
prob <- matrix(0,nrow(ci),nrow(cp))
for(k in 1:nrow(cp))
{
 for(i in 1:nrow(ci))
 {
  for(j in 1:ncol(ci))
  {
   if(datalist[i,1,k]>0 && datalist[i,2,k]>0)
   {
    index[i,k] <- 1
    prob[i,k] <- dbinom(i-1,n,cp[k,1])
   }
  }
 }
}

icr <- index*prob
cr <- cbind(cp,apply(icr,2,sum))
scr <- cr[sort.list(cr[,2]), ]
mincr <- scr[1,2] 
parameter <- c(scr[1,1],scr[2,1])

print(list("confidence coefficient for Agresti-Coull interval"=mincr,
"corresponding p"=parameter))